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Abstract 


We  applied  the  Wigner  Distribution  Function,  a  distribution  function  of  time  and 
frequency  based  on  an  initial  function  of  either  of  those  variables,  to  a  series  of  time 
based  correlation  functions.  These  time  based  correlation  functions  were  the  result  of  a 
1 -dimensional  free  particle  wave  packet,  the  reactant  wave  function,  which  had 
propagated  through  a  quantum  potential  well  and  then  had  components  of  the  reactant 
wave  function  that  exited  the  opposite  side  of  the  well  auto-correlated  in  time  with  a 
stationary  1 -dimensional  free  particle  wave  packet,  the  product  wave  function.  This 
process  was  undertaking  in  order  to  generate  a  3-dimensional  depiction,  in  time  and 
frequency,  of  the  reactant  wave  functions  interaction  with  the  quantum  potential  well. 
Fortran  77  code  was  utilized  to  generated  the  time  propagation  of  the  reactant  wave 
function  by  means  of  the  Split  Operator  Method,  which  was  given  the  following  initial 
set  of  conditions;  x0  —  —  20(Bohr  radii),  k0  —  3(atomic  units),  and  8  —  l(Bohr  radius). 
A  series  of  potential  wells  with  variable  depths  were  implemented  into  the  code.  The 
code  then  computed  the  correlation  in  time  of  the  exiting  reactant  wave  function  with  a 
stationary  wave  function  before  applying  the  Wigner  Distribution  Function.  When 
Wigner  Distribution  Function  was  applied  to  the  time  correlation  function  many 
recognizable  features  on  the  potential  well  were  observed  from  the  3-deminsional  plot 
generated  including  transmission  resonance  energy  levels.  The  classical  time  of  arrival 
was  also  captured  by  the  Wigner  Distribution  Function.  As  a  useful  tool  the  Wigner 
Distribution  Function  provides  more  insight  into  the  quantum  interactions  of  chemical 
reactions  in  terms  of  time  and  frequency  than  traditional  spectrographic  analysis. 
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A  WIGNER  DISTRIBUTION  ANALYSIS 


OF  ONE  DIMENSIONAL  SCATTERING 


I.  Introduction 


1.1.  Motivation 

In  an  overall  sense,  gaining  an  appreciation  for  how  elements  of  a  system  change 
over  time  is  the  fundamental  result  of  most  scientific  observation.  In  chemistry,  one 
monitors  the  populations  of  reactant  species  and  the  resulting  product  species  of 
compounds  over  the  course  of  a  chemical  reaction.  In  physics,  one  observes  the 
characteristics  (mass,  velocity,  etc)  of  an  element  of  a  system  and  monitors  how  those 
characteristics  change  as  the  system  evolves  in  time.  Typically,  these  temporal 
observations  are  bookended  observations  of  the  reaction  or  interaction  being  monitored. 
That  is  to  say,  we  generally  obverse  things  before  the  event  in  question  and  state  the 
properties  known  of  the  system  as  the  initial  conditions.  Once  the  event  is  question  has 
occurred  we  again  review  the  properties  of  the  system  and  state  that  these  are  the 
outcomes.  Development  of  techniques  and  procedures  to  aid  in  an  effort  to  view  and 
analysis  the  continuous  evolution  of  a  system  would  prove  invaluable  in  the  furtherance 
of  scientific  research. 

Of  immediate  importance  to  the  Air  Force  Office  of  Scientific  Research  is  the 
development  of  tools  and  technique  that  would  aid  in  the  development  and  refinement  of 
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Diode  Pumped  Alkaline  Lasers  (DPALs).  In  such  a  system,  an  Alkaline  atom  has  had  its 
electrons  pumped  into  the  P3/2  state  and  then  by  means  of  scattering  allowed  to  decay  to 
the  P 1/2  state.  Due  to  the  population  inversion  created  between  the  Pl/2  state  and  the 
ground  state,  a  situation  for  lasing  is  generated.  Fully  understanding  how  the  scattering 
works  in  the  decaying  from  the  upper  state  to  the  lower  one  may  aid  in  advancing 
techniques  for  making  such  Diode  Pumped  Lasers  more  efficient  and  feasible.  As 
mentioned  beforehand,  the  ability  to  monitor  such  interaction  over  time  in  a  useful 
manner  would  be  of  great  significance,  but  such  techniques  are  difficult  to  apply  due  to 
the  complexity  of  the  interactions  in  question  since  they  can  occur  in  multiple  dimensions 
and  amongst  several  dependant  variables. 

One  such  way  to  monitor  the  development  of  a  system  in  quantum  mechanics  and 
in  physical  chemistry  is  to  develop  a  correlation  function  in  time  of  two  wave  functions, 
one  the  reactant  wave  function  and  the  other  the  product  wave  function.  In  this  manner, 
the  reactant  wave  function  is  allowed  to  react  with  the  system  in  question  over  time  while 
being  correlated  with  a  stationary  product  wave  function  which  has  no  interaction  with 
the  system.  This  technique  has  many  benefits  as  it  records  the  information  of  the  reactant 
wave  function  over  time  and  eliminates  the  issues  imposed  by  multidimensional  or  many 
valued  problems  as  they  reduce  all  the  information  of  the  system  down  to  a  scalar 
function  of  time,  a  signal.  Once  a  signal  of  the  interaction  has  been  recorded,  it  is 
possible  to  apply  many  signal-processing  techniques  upon  it  to  regain  information  about 
the  initial  interaction. 
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One  such  signal  processing  technique  is  the  application  of  the  Wigner 
Distribution  Function  (WDF).  [1]  In  1932,  Eugene  Wigner,  in  an  attempt  to  return 
quantum  mechanics  to  the  domain  of  phase  space,  was  able  to  generate  a  function  of  both 
position  and  momentum  in  phase  space  from  the  wave  function  in  position  space  of  a 
quantum  state.  This  distribution1 2  function’s  integrated  projections  onto  either  the 
position  or  momentum  axes  returned  the  probability  densities  of  the  wave  equations  for 
those  variables  respectively.  The  WDF  is  given  by  the  form:  [2] 

1  f  2  ip y 

W(x,p)  =  —  \  xp*(x  +  y)xp(x  —  y)e  &  dy  (1-1) 

7rh  J 

where  x  is  position,  p  is  momentum,  and  xfi  is  the  time  independent  wave  function.  This 
phase  space  picture  of  quantum  mechanics  is  still  very  different  from  the  classical  phase 
space  picture  since  the  WDF  can  take  on  negative  values  and  is  a  surface  in  a  third  axis 
rather  than  a  parametric  plot  in  the  plane  of  the  position  and  momentum.  This  feature 
predominantly  arises  from  the  Heisenberg  Uncertainty  Principle,  which  makes  it 
impossible  to  simultaneously  know  both  the  position,  and  the  simultaneous  momentum, 
of  a  quantum  mechanical  system.  Figure  1  shows  a  WDF  of  the  initial  wave  function, 
ip(x,  0),  of  a  free  particle  wave  function  with  an  initial  position  of  0,  an  initial 
momentum  of  1,  and  an  initial  spread  of  1  . 


1  This  can  easily  be  confused  as  meaning  a  probability  distribution  function,  which  the  WDF  is  not  due  to 
the  fact  that  it  can  take  on  negative  values. 

2  The  meaning  behind  these  initial  conditions  and  the  equations  that  govern  them  will  be  illuminated  in 
Section  2.3,  but  for  completeness  that  are  included  here. 
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Figure  1.  The  WDF  of  the  initial  wave  function,  \p{x,  0),  of  a  free  particle  with  an  initial  position  of  0,  an  initial 
momentum  of  1,  and  an  initial  spread  of  1.  We  can  see  the  overall  Gaussian  feature  displayed  in  phase  space 
shifted  off  center  from  the  origin  in  the  momentum  axes.  This  is  due  to  the  wave  function  having  an  initial 
momentum  of  1  which  aligns  with  where  the  peak  of  the  Gaussian  is  located  in  phase  space  with  respect  to  the 
momentum  axis 


The  WDF  does  not  provide  any  new  information  than  what  was  already 
determined  by  the  individual  wave  equations.  However,  as  much  as  a  graph  does  not  tell 
you  anything  more  than  the  table  of  data  it  is  representing,  the  WDF  provides  a  new  way 
of  viewing  the  information  obtained  from  the  Schrodinger  Equation.  A  similar 
relationship  between  time  and  frequency  is  present  in  all  wave  mechanics,  as  is  displayed 
in  position  and  momentum  in  quantum  mechanics,  so  it  is  possible  to  extend  the  WDF  to 
the  time/frequency  domain.  Ultimately,  the  resulting  depiction  of  the  time  and  frequency 
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WDF  could  provide  insightful  representations  of  quantum  mechanical  events,  such  as  the 
scattering  of  a  wave  packet  upon  potential  wells  and  a  means  to  characterize  those  events. 

1.2.  Problem  Statement 

We  want  to  propagate  a  free  particle  reactant  wave  packet  through  a  quantum 
potential  well  and  compute  its  time  correlation  function  with  a  stationary  product  wave 
packet.  Then  we  want  to  apply  a  WDF  to  the  time  correlation  function  and  analyze  the 
results  using  known  attributes  of  the  quantum  square  well  that  the  reactant  wave  packet 
had  propagated  through.  Ultimately,  we  want  to  be  able  to  compare  these  results  with 
other  signal  processing  techniques,  such  as  a  spectrogram  of  the  same  signal,  and 
determine  if  the  WDF  is  a  viable  or  superior  method  of  analyzing  the  time  based  signal 
results  of  a  chemical  or  quantum  reaction.  Additionally,  we  want  to  determine  whether 
the  WDF  gives  us  any  further  or  new  insight  into  the  physics  within  the  well  past  the 
initial  scattering  of  the  reactant  wave  function. 

1.3.  Approach 

Solutions  to  the  equations  and  procedures  described  in  our  problem  will  be 
difficult  to  achieve  by  analytical  means  alone.  We  will  instead  employ  a  wide  variety  of 
numerical  techniques  in  order  to  simulate  our  wave  packet  scattering  event  in  a  computer 
model.  By  ensuring  that  we  take  sufficiently  small  numerical  increments  for  each  step, 
we  will  be  able  to  provide  the  best  approximation  of  the  true  analytical^  result. 

3  By  stating  “true  analytical  result’’  it  is  tempting  to  believe  I  am  referring  to  an  actual  analytical  answer  to 
our  problem  which  is  not  the  intent,  but  instead  to  highlight  the  analytical  result  that  would  have  been 
obtained  had  it  been  possible  to  take  the  considerable  amount  of  time  necessary  to  mathematically  construct 
a  solution  to  the  Hamiltonians  being  analyzed  here.  A  process  nature  accomplished  in  considerably  less 
time  than  men. 
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As  the  framework  for  our  computer  model  we  will  define  a  set  of  position 
coordinates  progressing  from  a  defined  —xmax  to  xmax  in  increments  of  Ax.  Using  this 
x-grid  we  will  discretely  define  a  potential,  V  (x),  such  that  for  all  values  within  a  and  -a 
of  the  origin  the  potential  is  a  set  negative  value  and  zero  for  all  other  locations.  This 
V (x)  will  serve  as  out  potential  well.  In  this  space  we  will  position  two  wave  packets  on 
either  side  of  the  potential  well  ensuring  that  the  wave  functions  are  a  sufficient  distance 
from  the  well  that  the  functional  values  will  have  sufficiently  decreased  to  zero  before 
occupying  the  same  space  as  the  well.  The  wave  function  to  the  left4  of  the  well  will 
serve  as  the  reactant  wave  function,  (x,  0),  for  our  numerical  wave  function 
propagator,  (x,  t),  while  the  wave  function  to  the  right  will  serve  as  our  product  wave 
function,  Tp(x,  0).  Figure  2  depicts  this  initial  setup. 


Figure  2.  The  initial  setup  of  the  reactant  wave  function  and  the  product  wave  function  positioned  opposite  each 
other  from  a  square  well.  In  this  depiction  the  wave  functions  are  simple  Gaussian  curves  centered  about  4 
(dashed  curve)  and  -4  (dotted  curve)  with  a  square  well  between  them  centered  about  the  origin  (solid  line).  The 
Gaussian  to  the  left  of  the  potential  well  is  the  reactant  wave  function  and  the  Gaussian  to  the  right  of  the 
potential  well  is  the  product  wave  function. 


4  By  left  we  mean  the  wave  function  is  centered  at  a  value  of  x  less  than  the  lower  bound  of  the  potential 
well.  The  inherent  symmetry  imposed  by  this  setup  is  not  an  imposed  restriction  of  our  numerical  method 
or  has  any  specific  benefit,  but  is  purely  a  matter  of  preference  and  the  methods  applied  could  be  applied  to 
any  1  dimensional  setup. 
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By  taking  the  correlation  of  these  two  functions  as  time  progress  as  the  reactant 
wave  function  propagates  towards  the  product  wave  function  on  the  other  side  of  the 
well,  we  will  record  the  amplitude  and  frequency  of  the  reactant  wave  function  after  it 
leaves  its  interaction  with  the  well.  This  is  possible  since  the  product  wave  function  is 
only  known  zero,  within  a  numerically  reasonable  degree5,  in  the  region  just  after  the 
well  in  relationship  to  the  reactant  wave  function  so  that  when  it  is  correlated  it  only  has 
non-zero  results  when  elements  of  the  reactant  wave  function  have  passed  through  the 
potential  well  and  are  then  propagating  through  the  region  where  the  product  wave 
function  is  located. 

Having  setup  our  initial  conditions,  we  will  then  develop  a  computer  algorithm  to 
replicate  the  Split  Operator  Method  of  propagating  the  initial  wave  function,  (x,  0), 
over  a  region  of  space  with  the  predefined  potential,  V  ( x ).  The  result  of  this  algorithm 
will  serve  as  the  propagated  wave  function  in  space  and  as  time  elapse  scattering  of  the 
wave  function  as  it  passes  the  potential  well,  V (x).  We  will  allow  the  propagator  to 
evolve  in  time  and  at  each  time  step,  we  will  generate  a  time  based  correlated 
function,  C  (t),  between  the  reactant  wave  function  and  the  product  wave  function,  by 
taking  the  integrated  product  of  the  two  over  position.  Equation  (1.2)  illustrates  this. 

C(t)  =  j  ipR*(x,t)ifjP(x,0)dx  (1.2) 

5  Technically,  the  function  is  never  truly  non-zero  as  it  is  continuous,  but  for  the  purposes  of  the  numerical 
approach  the  value  held  by  the  function  is  sufficiently  small  in  comparison  to  the  noise  generated  buy 
numerical  methods.  In  most  of  our  cases  the  functional  value  was  well  below  10A-10  when  we  say  it  is 
zero. 
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The  correlated  function  will  serve  as  our  time  signal  of  the  scattering  event.  Then  a 
numerical  WDF  will  then  be  apply  to  the  results  from  the  time  correlation  in  the 
following  fashion. 

1  f  2  j  O)  T 

W(t,a))——\  C*(t  +  r)C(t  —  r)e  &  dr  (1.3) 

7rh  J 

Prior  to  generating  the  WDF  of  the  correlation  function,  two  outcomes  will  be 
generated  to  verify  the  accuracy  of  the  numerical  WDF  generated.  First,  the  correlation 
function  will  be  multiplied  by  its  complex  conjugate  and  plotted  versus  time.  Second,  the 
correlation  function  will  have  an  Inverse  Fourier  Transform  applied  to  it  to  generate  the 
correlation  function  in  frequency  space,  C(co),  which  in  turn  will  be  multiplied  by  its 
complex  conjugate  and  plotted  versus  frequency.  These  plots  will  later  be  compared  to 
the  integrated  projections  of  the  WDF  of  the  correlation  function  in  order  to  verify  that 
the  WDF  generated  is  accurate  for  the  correlation  function  provided. 

The  correlation  function,  C(t),  will  then  be  integrated  in  a  computer  algorithm  to 
generate  the  WDF,  W ( t ,  of),  and  in  turn  this  will  be  the  WDF  of  the  scattering  event 
between  the  propagating  wave  packet  and  the  potential  well.  The  WDF  will  then  be 
plotted  versus  time  and  frequency.  To  verify  the  accuracy  of  the  transform,  W (t,  co)  will 
then  be  integrated  versus  time  to  produce  the  WDF  projection  in  frequency  space,  and 
integrated  versus  frequency,  to  produce  the  projection  in  the  temporal  domain.  These 
projections  will  be  plotted  and  compared  to  the  plots  of  |C(co)|2  and  |  C  (t)  | 2  generated 
previously.  The  two  plots  should  appear  to  over  lap  for  the  WDF  to  be  accurate.  Figure 


3  shows  this  comparison  of  \xp(x)\2  and  \ip(k)\2  for  the  previously  shown  example  of  the 
WDF  in  figure  1  of  the  free  particle. 


Figure  3.  The  comparison  of  the  probability  distribution  in  both  position  and  momentum  of  free  particle  wave 
function  at  The  solid  red  lines  represents  the  solution  derived  by  taking  the  modulus  of  the  wave 

function  and  the  blue  dots  represent  the  solution  of  the  integrate  projection  of  the  WDF  with  the  axes  of  interest. 

With  this  approach  we  will  first  analyze  the  WDF  for  a  square  potential  well.  By 
using  the  square  well  we  can  compare  the  known  attributes  with  the  results  of  the  WDF 
and  determine  if  the  WDF  provides  any  insight  into  the  physics  going  on  within  the 
scattering  event.  Additionally,  we  can  compare  the  WDF  result  with  a  spectrogram  of  the 
correlation  function. 

In  Section  II,  we  will  discuss  giving  a  broader  scope  of  the  mathematical 
techniques  used  to  generate  the  WDF  such  as  the  Fourier  Transform,  Split  Operator 
Method,  and  the  Free  Particle  Wave  Function.  A  treatment  of  the  square  well  will  also  be 
included  for  the  purpose  of  comparison  to  the  WDF.  Section  HI  will  go  over  the 
technical  aspects  of  the  computer  algorithm  generated  for  this  research  in  order  to  give  an 
understanding  of  how  the  techniques  were  used  and  an  explanation  of  the  parameter 
decisions  made.  Section  IV  presents  the  results  of  the  generated  WDF  for  the  various 
wells  analyzed,  followed  by  an  analysis  in  Section  V. 
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II.  Background 


2.1.  Fourier  Transform 

A  crucial  component  of  research  between  position  and  momentum  in  quantum 
mechanics,  as  well  as  time  and  frequency,  is  the  proper  use  of  Fourier  Transformations. 
Fourier  Transformation  is  a  process  which  transforms  a  complex-valued  function  of  a  real 
variable  into  another  variable.  [3]  [4]  A  better  way  to  look  at  this  transformation  is  by 
viewing  a  function  of  any  variable  as  vector  of  some  magnitude  and  direction  in  an 
infinitely  defined  space  of  the  variable.  So  if  f{xf)  is  the  value  of  the  function  at  x1  then 
in  this  vector  view,  f(xf)  is  the  projection  of  /  on  the  xt  axis.  Extending  this  logic  to  a 
continuous  function  which  existed  over  the  entire  x  domain,  we  would  be  able  to  describe 
a  vector/ to  be  in  the  infinite  vector  space  of  x.  If  this  vector  was  then  transported  in  a 
new  vector  space,  with  each  new  axis  being  related  to  a  corresponding  x-space  axis  by  a 
phase  shift,  the/vector  would  have  a  new  set  of  projections  onto  the  new  infinite  vector 
space  and  would  have  a  new  shape  in  the  new  domain. 

This  is  essentially  how  the  Fourier  Transform  works  between  conjugate  pairs  of 
variables.  For  the  conjugate  pair  of  position  and  momentum;  its  relationship  is  given  by 
the  following:  [5] 

f{x)  =  —=  [  F(k)elkxdk  F(k)  =  —=  f  f(x)e~lkxdx  (2.1) 

y/2n  J  yJ2nJ 

1 

The  factor  ~^=  is  used  to  normalize  the  functions  as  you  transform  between  one  domain, 

to  the  other,  and  back  again.  With  that  in  place  you  can  go  back  and  forth  between  the 
two  domains  and  always  return  to  the  original  functions. 
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Another  conjugate  pair  is  time  and  angular  frequency  which  is  given  by  the  same 
form  as  the  equation  (2.1)  with  the  substitution  of  t  for  x  and  or  for  k.  By  extension,  this 
time/angular  frequency  pairing  is  also  a  time/energy  pairing  due  to  the  DeBroglie 
relationship,  E  —  ho>. 

In  order  to  analytically  compute  the  Fourier  Transform,  the  function  must 
converge  as  the  domain  goes  to  its  infinite  limits.  This  poses  a  problem  when  you  look  at 
things  in  a  discrete  and  finite  frame  of  view.  However,  as  long  as  the  function  value 
sufficiently  goes  to  zero  well  before  the  end  of  the  finite  domain,  a  discrete  form  of 
equation  will  provide  the  correct  transform  for  the  original  function. 

2.2.  Wigner  Distribution  Function 

The  development  of  quantum  mechanics  and  its  counter  intuitive  nature  posed 
many  problems  for  physicists  in  the  early  part  of  the  1900s.  Classically,  a  system  would 
be  modeled  and  evaluated  based  on  known  characteristic  values,  such  as  the  precise 
position  of  a  billiard  ball  at  an  exact  time,  and  from  that  information  future  outcomes 
could  be  stated  with  a  fair  amount  of  certainty  based  on  the  laws  of  classical  dynamics.  If 
you  knew  the  position  and  momentum  of  a  particle  and  the  system6  it  occupied  then  you 
knew  not  only  everything  that  had  happened  to  the  particle,  but  what  would  happen  to  it 
given  a  universal  knowledge  of  the  system  it  occupied.  Laplace  is  famously  noted  as 
having  stated  that  the  universe  was  composed  of  particles  obeying  Newton’s  laws  and 
“once  the  interaction  between  these  bodies  is  precisely  known  and  the  position  and 
velocities  of  all  the  bodies  at  any  given  instant  are  known,  these  coordinates  and 

6  This  state  meant  makes  the  broad  assumption  that  the  system  is  a  non-interactive  closed  system  with 
regard  to  the  rest  of  the  universe.  More  of  a  mental  exercise  then  an  actual  case. 
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velocities  are  determined  (through  Newton’s  second  law)  for  all  time.”  [6]  This 
comfortable  determinism  would  soon  be  shaken  by  DeBroglie’s  Wave-particle  Duality 
and  Heisenberg’s  Uncertainty  Principle. 

However,  this  did  not  discourage  some  from  approaching  the  new  landscape  of 

7 

physics  with  some  traditional  notions.  One  of  these  physics  was  Eugene  Paul  Wigner 
who  tried  to  tie  in  the  traditional  phase- space  picture  with  what  was  being  developed  in 
the  field  of  quantum  mechanics  and  its  force  condition  of  quantum  uncertainty.  In  1932, 
he  developed  a  distribution  function  as  a  means  to  study  quantum  corrections  to  classical 
statistical  mechanics.  Fundamentally,  this  quasi-distribution  function  generates  a  2D 
surface  in  phase-space  whose  integrated  projection  on  either  axis  will  return  the  modulus 
of  the  wave  function  for  that  axis.  Today  this  distribution  function  is  called  the  Wigner 

o 

Distribution,  or  the  Wigner-Ville  Distribution  ,  and  is  shown  as  equation  (2.2)  for  the 
conjugate  pairs  of  position  and  momentum. 

ivfep)  =TJ  W  +  yMx-y^dy  (2.2) 

Fundamentally,  the  distribution  function  is  a  Fourier  Transform  of  a  correlated  wave 
function.  An  example  of  such  a  distribution  function  is  given  in  figure  1.  In  figure  1,  all 
the  values  of  the  distribution  are  positive,  but  it  is  entirely  possible  for  the  WDF  to  render 
negative  values  in  the  plot.  The  meaning  behind  this  negative  property  of  the  distribution 


7  Pal  Jeno  Wigner,  his  native  Hungarian  name,  was  born  on  November  17,  1902  in  Budapest.  He 
immigrated  to  the  United  States  in  1930.  He  was  a  heavy  contributor  to  Quantum  Group  Theory.  He  died 
on  January  1,  1995.  His  notable  relative  is  Paul  Dirac,  his  brother-in-law  by  his  first  marriage. 

8  After  J.  Ville  who  in  1948  rederived  it  as  a  quadratic  (in  signal)  representation  of  the  local  time -frequency 
energy  of  a  signal.  It  is  additionally  interesting  to  note  that  both  Heisenberg  and  Dirac  also  derived  this 
function,  but  missed  its  significance.  [21]  Also  in  1949,  Jose  Enrique  Moyal  independently  derived  it.  [20] 
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is  unclear  and  makes  it  exceedingly  difficult  for  one  to  interpret  the  WDF  as  a  probability 
distribution  in  phase  space. 

Additionally,  figure  2  depicts  the  comparison  of  the  integration  of  the  WDF  along 
the  momentum  axis  with  the  modulus  wave  function  in  its  coordinate  representation. 

This  demonstrates  that  the  WDF  contains  the  information  about  the  wave  function  in  both 
coordinate  space  and  momentum  space. 

So  far  we  have  examined  the  position-momentum  WDF,  but  the  WDF  is  not 
limited  to  that  singular  case.  Any  conjugate  pair  can  have  this  technique  applied  to  it  and 
it  is  often  used  in  the  pairing  of  time  and  frequency  such  as  shown  in  equation  (2.3). 

1  f  2  i  a)  t 

W(t,a))=—\  ip*(t  +  z)xp(t  —  z)e  ft  dz  (2.3) 

7rh  J 

This  is  the  form  of  the  WDF  we  will  use  when  examining  the  correlation  of  out  scattered 
wave  function  over  time.  Except  in  our  case  ip(t)  in  equation  (2.3)  will  be  replaced  by 
the  calculated  time  correlation  of  the  reactant  and  product  wave  function,  C(t),  discussed 
later  in  Section  2.6. 

2.3.  Propagating  Wave  packet 

Foremost  in  our  analysis  of  the  WDF  as  a  signal  processing  tool  is  ability  to 
accurately  propagate  our  reactant  wave  function.  As  this  function  will  serve  as  our 
propagating  element  in  our  initial  setup  and  will  be  the  wave  function  that  will  scatter 
with  the  quantum  potential  well  in  order  to  later  produce  the  time  signal  by  means  of 
correlation  with  our  product  wave  function.  The  uniqueness  in  the  expression  of  the  free 
particle  wave  packet  is  that  as  a  wave  function  the  function  needs  to  be  continuous  over 
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all  space,  but  still  maintain  its  relative  localization  as  it  propagates  in  time  much  like  a 
particle.  Many  authors  discuss  the  development  of  this  expression.  [7]  [6] 

Schrodinger’s  Equation  tells  us  that  the  time  evolution  of  a  wave  function  is  the 
product  of  the  initial  wave  function,  an  eigenfunction  of  the  Hamiltonian,  with  the 
exponential  of  the  Hamiltonian  over  time. 

(2.' 

-i  H  t 

ip(x,t)  =  e  fr  ip(x,  0) 

Recall  from  our  discussion  on  Fourier  Transforms  that  a  function  of  position  can  be 
defined  as  the  integration  of  a  function  in  momentum  space,  as  shown  in  equation  (2.1), 
and  then  equation  (2.4)  can  be  rewritten  as 


-iHt  1  r  (2-5) 

t)  —  e  &  —=  I  F(k)elkxdk 
\]2n  J 

Since  the  free  particle  assumes  that  the  potential  component  of  the  Hamiltonian  is  zero 
over  all  space,  only  the  kinetic  part  remains  which  is  also  a  function  of  momentum, 

(hk)2 

T  —  =  h co.  Equation  (2.5)  then  consolidates  into  the  following 


1  f  (2-6) 

ip(x,t)  ——=  I  F(/c)el(fcx  wt^dk 
V27T  J 

However,  F(k),  is  time  independent  and  by  the  same  Fourier  relationship  can  be  given  by 
the  following 


F(k)  = 


(2.7) 
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Substituting  equation  (2.7)  into  equation  (2.6)  we  reach  the  following  expression 

<k*,o  =  2- 

Since  we  want  our  initial  wave  function  to  have  a  relatively  localized  position  we 
will  use  the  time  independent  Gaussian  wave  function 


"  r 

j  * 


ip(x',  0)e  lkx' dx'  el(kx  cot^dk 


ip(x,  0)  = 


S1/2(2ny /4 


eik0xe~x2/4^ 


(2.9) 


as  our  initial  wave  function  in  equation  (2.8)  and  reach  the  following  form 


t)  =  _i_  [  e('^W  f 


(2.10) 


<52(271)4 


where  t  =  .  We  now  integrate  equation  (2.10)  over  all  space  and  all  momentum  to 

arrive  at  out  final  form  of  the  free  particle  wave  function 


ip(x,t)  = 


(x-x0)2  .  .  2 

4  §2  +  1  Ko  ^  2  mK° 


1  + 


\l2n  S  2 


i  h  t 
2m62 


1  + 


i  h  t 
2  m  82 


(2.11) 


where  x0  is  the  initial  center  position  of  the  wave  packet,  k0  is  the  initial  kinetic  energy, 


8  is  the  initial  spread  of  the  wave  packet  in  coordinate  space,  and  m  is  the  mass  of  the 
particle.  Figure  4  depicts  the  time  evolution  of  the  real  and  imaginary  component  of  the 
wave  function. 
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t  =  1  (atomic  time  units)  t  =  5 (atomic  time  units) 


Figure  4.  The  time  evolution  of  the  real  and  imaginary  components  of  a  free  particle  with  initial  values  of  x  =  - 
20,  k  =  3,  and  a  spread  of  1. 

As  you  can  see,  the  wave  function  starts  with  a  distinct  position9  and  spreads  out 
over  all  space  as  time  progresses.  This  is  the  expected  outcome  of  the  time  evolution  of 
the  free  particle.  Initially,  the  probability  of  finding  the  particle  in  a  set  region  is  high, 
but  as  time  progress  the  probability  of  location  of  the  particle  gets  spread  over  an  ever 
larger  location.10  When  we  set  t  =  0  equation  (2.11)  simplifies  down  into, 


ip(x,  0)  = 


Vl  n  82 


(x-x0)2 

4S2 


+  i k0 ( x-x0 ) 


(2.12) 


9  Saying,  “distinct  position,”  may  be  to  strong  of  a  statement.  The  particle  location  is  as  distinct  as 
Heisenberg’s  Uncertainty  principle  allows,  but  in  terms  of  the  entirety  of  space,  the  location  is  very 
localized. 

10  Actually,  the  particle  has  a  probability  of  being  located  at  any  point  in  the  continuum  of  position  space, 
but  the  beauty  of  the  free  particle  wave  function  is  that  the  overall  Gaussian  shape  restricts  the  region  of 
highest  probability  over  a  smaller  region  which  only  broadens  with  time. 
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Equation  (2.12)  is  exactly  equation  (2.9)  which  you  may  recall  was  the  time  independent 
Gaussian  wave  function.  This  is  not  surprising  as  when  you  take  the  time  dependent 
expression  of  the  wave  function  and  set  t  =  0  you  should  return  to  the  time  independent 
form  of  the  expression.  This  form  will  serve  as  our  initial  wave  function,  but  the  time- 
dependent  form  will  not  function  as  our  time  propagator  because  it’s  Hamiltonian  does 
not  contain  the  elements  of  the  potential  well  we  wish  to  have  the  particle  scatter  upon  as 
explored  above. 

2.4.  Split  Operator  Method 

Previously  we  demonstrated  the  form  of  the  propagating  wave  packet  wave- 
function  when  the  Hamiltonian  has  no  potential  energy  component.  However,  in  our  case 
we  will  be  looking  at  examples  where  the  potential  is  not  zero  everywhere.  This  would 
force  us  to  re-exam  our  previous  derivation  for  each  Hamiltonian  we  wish  to  look  at  in 
order  to  derive  the  appropriate  eigenfunction  that  showed  the  correct  time  based 
propagation.  As  we  wish  to  exam  a  wide  variety  of  wells,  this  would  become  a  very 
laborious  process  and  with  possibility  of  finding  no  analytic  solution. 

Instead  of  developing  a  correct  function  for  each  Hamiltonian,  we  will  be  using 
the  Split  Operator  Method  to  develop  a  time  based  wave  function  based  on  an  input 
potential.  The  Split  Operator  Method  [8]  [9]  is  an  outcrop  of  the  Schrodinger  Equation  as 
show  in  equation  (2.5).  Split  Operator  Method  takes  the  time  dependent  form  of  the 
Schrodinger  Equation  and  rewrites  the  Hamiltonian  Operator  in  terms  of  its  component 
operators,  H  =  T  +  V.  In  this  new  form  we  can  rewrite  the  argument  of  the  exponential 
and  a  new  form  of  the  Schrodinger  Equation  appears  as:  [10] 
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(2.13) 


i  V  At  -i  T  At  i  V  At 

ip(x,t  +  At)  =  e  2h  e  h  e  2h  ifj(x,  t)  +  0(At3) 

The  factor  of  0(At3)  arises  from  the  fact  that  T  and  V  do  not  commute  and  the  third  term 
in  equation  (2.14)  is  not  1  as  it  is  when  applied  to  scalars.  [11] 


gCib  _  qQ. gb q  (2.14) 

In  order  to  avoid  issues  with  this  0(At3)  term,  we  will  ensure  that  we  make  At 
sufficiently  small  in  comparison  to  the  total  energy  of  the  system  we  want  to  examine  and 
we  will  assume  that  the  last  term  is  small  enough,  by  the  third  power,  to  be  considered 
zero.  With  the  Schrodinger  Equation  now  in  the  form  given  by  equation  (2.13),  we  can 
capitalize  on  the  Fourier  Transform  to  avoid  further  complexities  of  the  T  and  V 
operators.  By  this  I  mean  in  coordinate  space,  the  V  operator  is  merely,  V ( x )  or 
diagonalized  in  coordinate  space,  which  makes  the  operation  of 


t  V(x)  At 

f  (x,  t)  =  e  2h  ip(x,t) 


(2.15) 


a  product  of  two  functions.  However,  in  coordinate  space  T  becomes  harder  to  manage 
as  shown  in  equation  (2.16) 


■  s2  . . 
ls^At 


(2.16) 


h(x,  t)  —  e  h  f(x,  t) 

To  avoid  expanding  the  exponential  into  a  Taylor  Series  and  then  applying  the  operator  to 
the  function,  we  can  simply  Inverse  Fourier  Transform,  f(x,  t),  into  momentum  space 
and  make  the  result  a  simple  product  of  two  functions  as  before.  As  T  in  momentum 
space  is  simply, 
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(2.17) 


T  =  T(k ) 


(hfc); 

2m 


which  is  to  say  T  is  diagonalized  in  momentum  space  and  then  equation  (2.16)  can  be 


rewritten  as: 


i  T(k)  At 

h(k,  t)  =  e  h  /(fc,  t) 


(2.18) 


From  that  point  we  simply  Fourier  Transform  back  to  the  coordinate  representation  and 
compute  the  last  product  to  arrive  at  xp(x,  t  +  At).  In  repeating  this  process  by  using  the 
resultant  i//(x,  t  +  At)  as  the  new  ip(x,  t),  the  time  evolution  of  the  wave  function  can  be 
generated  in  steps  of  At.  Equation  (2.19)  shows  the  overall  process: 


-i  V(x )  At 

xp{x,  t  +  At)  =  e  2l  FT 


-i  T(k)  At 

■  -iV(x)At 

e  h  FT-1 

e  2h  xjj(x,t) 

(2.19) 


2.5.  Square  Well  Scattering 

Understanding  of  the  quantum  square  well  is  crucial  to  our  analysis  of  the  WDF 
of  a  scattered  wave  packet  incident  upon  the  well.  Fortunately,  the  finite  square  well  is  a 
well  know  problem  with  multiple  examples  to  pull  from.  [12]  [6]  [13]  [7]  By 
understanding  its  known  attributes  we  will  be  better  able  to  identify  features  of  our 
resultant  WDFs  and  make  knowledgeable  observations. 

If  we  examine  the  scenario  where  we  have  a  propagating  wave  which  travels  in 
the  direction  of  x  and  experiences  zero  potential  until  a  region  of  space,  from  -a  to  a, 
where  the  potential  is  a  negative  constant  value,  we  then  have  the  case  of  a  quantum 
square  well.  If  we  accept  that  a  portion  of  the  propagating  wave  will  reflect  off  the  first 
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boundary  well  and  that  some  will  pass  through,  we  can  then  define  the  equation  of  the 
wave  in  the  region  before  the  well  as: 

x />(*)  =  Aeikx  +  Be~ikx  (2.20) 

V2  m  E  (2.21) 

where  A  is  the  incident  amplitude  and  B  is  the  reflected  amplitude.  Within  the  well  we 
can  define  the  wave  function  as: 


—  C  sin  lx  +  D  cos  lx 

j2m(E  +  V0) 
h 

On  the  far  side  of  the  well  we  can  define  the  wave  function  as: 


(2.22) 

(2.23) 


xjj(x)  =  Felkx  (2.24) 

In  this  case,  k  is  the  same  as  equation  (2.21)  and  there  is  no  reflected  amplitude  as  we 
assume  the  wave  is  travelling  from  left  to  right.  With  these  equations  we  impose  two 

boundary  conditions  on  the  wave  function.  First  that  xp(x )  is  continuous  and  that  ^  also 


is  continuous  at  the  boundaries. 


From  that  we  can  set  up  a  system  of  equations  that  allows  us  to  set  the  transmitted 
amplitude,  F,  in  terms  of  the  incident  amplitude,  A. 


-2  ika 


F  = 


cos  2 la  —  i  ^2fc/a  +  l2) 


A 


(2.25) 
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We  can  then  define  the  transmission  coefficient  of  the  well  by  using  equation  (2.25)  in 


equation  (2.26). 


(2.26) 


2.5.1.  Transmission 

The  transmission  coefficient  is  a  function  of  energy;  when  equation  (2.25)  is 
applied  to  equation  (2.26)  what  is  left  is  an  expression  in  terms  of  k  and  1,  which  they 
themselves,  as  show  in  equation  (2.21)  and  (2.23),  are  functions  of  energy.  Equation 
(2.27)  demonstrates  this. 

T  = _ 1 _  (2-27) 

1  +  4huj(hL  +  Va) sin2  (lT ^2m<ha>  +  V) 

In  figure  5,  we  see  this  transmission  coefficient  plot  over  energies  of  the  incident  wave 
function  and  we  can  see  that  when  energy  increases,  the  degree  of  the  transmission  is 
increased. 
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Figure  5.  The  transmission  as  a  function  of  energy  of  a  square  well  with  depth,  V„  =  20,  and  width,  2 a  =  20. 


You  will  also  notice  the  sinusoidal  shape  of  the  curve  as  it  moves  towards  1.  These  peaks 
rise  up  to  the  value  of  total  transmission  before  returning  to  the  general  arc  of  the  curve. 
Since  E  is  related  to  angular  frequency  by  h,  these  points  represent  frequencies  of  total 
transmission  through  the  square  well.  At  these  frequencies/energies  the  square  well  is 
completely  transparent  to  the  particle. 

The  overall  shape  of  the  transmission  coefficient  is  the  same  for  all  square  wells, 
but  varies  in  the  degree  at  which  the  peaks  rise  up  from  the  general  arc  of  the  function  to 
1  and  in  the  spacing  between  these  peaks  (i.e.  the  spacing  between  peaks  decreases  as  the 
well  depth  increases).  In  Merzbacher’s  analysis  [7]  of  the  square  well,  he  defines  a  /? 
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parameter  to  the  square  well  that  determines  the  degree  of  resonance  given  by  the 
following: 


'2mVn 


(2.28) 


More  resonate  wells  are  seen  to  have  a  larger  /?  value  and  thus  has  a  larger  depth  or 


width,  or  the  particle  has  a  larger  mass,  or  a  combination  of  all  three. 


2.5.2.  Energy  Levels 

Additionally,  the  peaks  of  the  transmission  coefficient  coincide  with  energies  that 
cause  the  sine  squared  term  to  be  a  multiple  of  n. 


(nn)2h  V0  (2.30) 

2  m  (2  a)2  +  h 


2.6.  Time  Correlation  Function 

So  far  we  have  examined  the  propagation  of  a  reactant  wave  function  in 
conjunction  and  its  scattering  off  a  potential  well,  but  not  a  means  of  recording  the  result 
of  the  scattering  as  a  function  of  time.  For  that  we  must  correlate  our  wave  function  with 
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a  stationary  wave  function  at  a  given  position  beyond  the  well  from  the  starting  point  of 
our  propagating  wave  function.  [14] 

C(t)  =  J  ipR*(x,t)iftP(x,0)dx  (2-31) 

The  correlation  function  enables  us  to  see  what  elements,  the  amplitudes  of  the  real  and 
imaginary  components,  of  the  propagating  wave  function  are  passing  a  certain  position  as 
time  evolves.  Additionally,  we  will  also  gain  information  about  frequencies  that  pass 
through  the  stationary  point  based  on  the  Fourier  relationship  between  time  and 
frequency. 

2.7.  Atomic  Units 

As  we  explore  the  interactions  of  the  propagating  wave  function  with  the  imposed 
well  we  will  use  the  atomic  unit  system  for  our  numerical  values.  This  is  a  common  unit 
system  for  analysis  of  events  at  the  atomic  level  and  is  often  used  in  Chemical  Physics 
research.  The  main  benefit  it  will  serve  in  this  research  is  in  reducing  many  of  our 
constants  values  to  1  such  as  mass,  the  mass  of  an  electron  in  this  case,  and  Planck’s 
constant.  While  the  numerical  approach  is  indifferent  to  the  units  used,  the  atomic  unit 
system  gives  us  a  more  intuitive  scale  upon  which  to  view  the  results  of  our 
computational  model.  Table  1  shows  a  comparison  of  the  units  in  the  atomic  unit  system 
with  those  in  the  International  System  of  Units  (SI) 
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Table  1.  Atomic  Units  System  [15] 


Quantity 

Name 

Symbol 

SI  Value 

mass 

electron  rest  mass 

mP 

9. 109*  10'31  kg 

length 

Bohr  radius 

h 

a0  — 

mPca 

5.291*10  11  kg 

charge 

elementary  charge 

e 

1.602*10‘iy  C 

angular  momentum 

Reduced  Planck’s 

constant 

h 

1.504*  10"34  J  s 

energy 

Hartree  energy 

Eh  —  mPc2a2 

4.359*10‘1S  J 

electrostatic  force 

constant 

Coulomb’s  constant 

1 

4tren 

8.987*  109  N  m2 

time 

h 

E~h 

2.418*10"17  s 

velocity 

&0  Eh. 

_ h _ 

2. 1 87*  106  m  s'1 
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III.  Numerical  Approach 


3.3.  Fast  Fourier  Transform 

Previously,  we  discussed  the  process  of  Fourier  Transforming  one  continuous 
function  in  a  position  (or  time)  domain  into  a  continuous  function  in  a  momentum  (or 
frequency/energy)  domain.  However,  since  we  are  taking  a  numerical  approach  to  the 
propagation  of  a  free  particle  scattering  off  of  a  square  well,  we  will  not  have  a 
continuous  domain  to  work  with,  but  the  inside  of  a  very  controlled  finite  domain.  This 
will  force  us  to  perform  a  Discrete  Fourier  Transform  when  we  seek  to  perform  the 
operations  given  by  equation  (2.19).  All  numerical  computation  where  programmed  into 
FORTRAN77. 


It  would  be  possible  to  take  the  integral  form  shown  in  equation  (2.1)  and  revert  it 
into  a  Riemann  Sum  expression  with  a  Ax  of  our  choosing 


"  IL 

F(kd  =  ^ 


Ax  *  /(x)  *  e 


iktx 


(3.1) 


X=X± 

This  would  be  able  to  discretely  transform  f(x)  on  the  range  x1  ->  xn  into  F(/q).  This 
would  be  repeated  for  every  F(/q)  over  a  pre-determined  range  and  would  provide  the 
Fourier  Transformed  function.  Provided  that  the  ranges  chosen  are  sufficiently  large 
such  that  non-zero  portions  of  the  functions  are  not  excluded  in  either  domains,  then  this 
Discrete  Fourier  Transform  will  be  highly  accurate. 


This  approach  however  has  a  drawback  when  it  comes  to  the  computational  time 
it  takes  to  perform  these  procedures.  “The  interpolation  of  2m  data  points  by  the  direct 
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2  9 

calculation  technique  requires  approximately  (2m)  multiplications  and  (2m)  additions.” 
[16]  Fortunately,  J.W.  Cooley  and  J.W.  Tukey  developed  a  numerical  process  that 
reduced  the  number  of  multiplications  and  additions  to  (m  log2  m),  called  the  Fast 
Fourier  Transform.  They  only  requirement  imposed  by  this  technique  is  that  the  position 
(time)  domain  and  the  momentum  (frequency)  domain  must  share  the  same  number  of 
grid  points,  and  those  grid  points  have  to  be  on  the  order  of  2N.  [16]  [3] 

Due  to  its  efficiency  there  are  many  available  Fast  Fourier  Transform  subroutines 
available  in  most  every  computational  programming  language.  We  implemented  a  1 
dimensional  Fast  Fourier  Transform  into  our  code  in  order  to  perform  the  transforms 
need  in  equation  (20)  and  also  for  additional  transforms  into  momentum  space  for 
calculating  the  norm  of  k  at  each  time  step  in  our  program. 

3.4.  Aliasing 

One  of  the  issues  encountered  with  our  use  of  discrete  FT  and  the  process 
efficient  FFT  was  one  of  Aliasing.  Aliasing  essentially  occurs  when  a  continuous  signal 
is  sampled  and  reconstructed  on  a  fixed  range  allowing  for  the  development  of  distortions 
or  artifacts  to  be  represented.  In  a  purely  mathematical  FT,  the  integration  of  the  input 
signal  would  be  over  the  entire  time  domain,  i.e.  from  t  =  —  oo  to  t  —  oo.  However, 
doing  such  integrations  numerically  is  impossible,  so  we  imposed  the  restriction  that  over 
a  short  window  the  signal  would  be  zero  at  the  imposed  boundaries  of  the  window  and 
therefore  the  integration  could  be  imposed  over  those  limits.  While  these  imposed 
restrictions  do  work  to  transform  the  function  from  time  to  frequency,  they  have  to  be 
handled  properly  or  they  will  create  distortions  or  artifacts. 
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The  reason  is  that  when  this  windowed  regime  of  the  function  is  generated,  it 
actually  creates  an  infinite  series  of  these  windows  that  follow  one  another  in  the  domain 
of  the  space  being  transformed  into.  Each  of  these  windows  is  a  replicate  of  the  previous, 
starting  where  the  previous  window  ended.  If  an  arbitrary  length  is  chosen  over  which  to 
measure  and  describe  the  domain  being  transferred  into,  then  there  is  the  risk  of  recording 
more  than  one  of  these  windows  and  creating  a  transform  that  actually  contains  multiples 
of  the  true  transform. 

These  distortions  can  be  avoided  by  imposing  a  range  restriction  on  the  domain 
being  transformed  into.  This  restriction  is  called  the  Nyquist  Frequency  and  is  shown  in 
equation  (3.2).  [4]  [3]  [17] 

n 

^ max  (3-2) 

This  restriction,  while  aiding  in  eliminating  distortions  and  artifacts,  forces  us  to  be  very 
careful  in  our  choice  of  At.  We  are  already  forced  to  choose  a  At  that  is  small  in 
comparison  to  the  energies  of  the  system  we  are  investigating  due  to  the  use  of  the  Split 
Operator  Method,  but  now  we  must  be  careful  not  to  make  our  0)max  too  large  or  the  A co 
too  big  to  be  a  meaningful  representation  of  \ip(a))\2,  but  also  ensuring  comax  included  all 
of  the  information  of  |i/;(co)  |2.  In  a  case  involving  x  ->  k  transforms,  we  are  further 
burdened  by  the  FFTs  requirement  that  both  domains  be  characterized  by  2N  number  of 
points  to  function  properly.  This  adds  the  further  balancing  act  of  ensuring  kmax 
includes  all  the  appropriate  information  of  the  kinetic  energy  from  the  scattering,  even  in 
addition  to  wanting  both  Ax  and  A k  to  be  small  enough  to  properly  characterize  the 
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functions  in  those  domains.  Table  2  shows  this  imposed  restriction  and  its  resulting 
values  used  in  the  FORTRAN  77  code  generated  for  this  research. 


Table  2.  Comparison  of  grid  sizes  and  grid  step  sizes  for  Discrete  Fourier  Transforms 


Conjugate 
Pairs 
(a  -»  b ) 

Ct-max 

A  a 

b 

umax 

Ah11 

#points 

x  ->  k 

100 

.09766 

32.16991 

.03142 

2048 

t  — >  (i) 

99.99 

.09 

34.90659 

.03142 

mi 

It  should  be  noted  that  At  was  actually  .01  in  the  Split  Operator  Method  portion  of 
the  code  but  was  down  sampled,  as  explained  previously,  to  .09  to  save  time  in  the 
processing  of  the  WDF  portion  of  the  code.  This  value  of  .09  served  as  the  value  of  At 
for  all  transforms  of  t  ->  co  in  the  code  generated  and  is  why  it  is  displayed  as  the  delta 
term  in  the  table  above. 

3.5.  Absorbing  Boundaries 

Since  we  are  performing  discrete  transforms  over  a  fixed  regime  of  space,  we  run 
a  similar  risk  to  the  aliasing  described  previously.  If  the  wave  function  is  allowed  to 
propagate  to  the  point  where  elements  of  the  function  pass  the  right  side  of  our  imposed 
window,  we  will  see  the  return  of  those  elements  at  left  side  which  will  interact  with 
elements  still  in  the  window  and  distort  our  wave  function.  In  order  to  avoid  this  we  need 
to  ensure  that  the  wave  function  goes  to  zero  as  it  reaches  the  end  of  the  special  regime. 
To  do  this  we  will  add  an  imaginary  term  to  our  potential  function  that  will  be  localized 
to  the  edges  of  our  special  window.  [18] 


11  It  was  not  intended  for  Ak  and  Aa>  to  match,  but  was  the  unexpected  outcome  of  the  values  used. 
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^total(X)  ^int C^O  "b  iVabs(.%) 


(3.3) 


When  this  term  is  placed  into  the  potential  exponential  of  equation  (2.19),  we  see  the 
following  result 


iVAt  iVrotaiAt 

e  2  —  e  2  ip(x,t ) 

(3.4) 

iV  At  i(V int+iV abs)^ 

2  l/j(x,  t)  —  e  2  xjj(x,  t) 

(3.5) 

iV  At  dVint~V  abs)^t 

2  lp(x,  t)  =  e  2  lp(x,  t) 

(3.6) 

iV  At  iV  in^At  V absAt 

2  xjj(x,  t)  —  e  2  e  2  x/j(x,t ) 

(3.7) 

Provided  that  Vabs(x)  is  negative,  the  second  term  of  equation  will  serve  to 
dampen  the  wave  function  multiplied  with  it.  Additionally,  to  ensure  the  absorbing 
boundary  potentials  do  not  drastically  dampen  the  wave  function  to  the  point  of  altering  it 
within  the  special  window,  it  is  best  to  use  a  function  that  gradually  applies  the  negative 
potential.  We  used  a  simple  Gaussian,  depicted  in  equation  (3.8),  to  accomplish  this. 

fX—Xmin\2  (X—Xmax\2  (3.8) 

Vabs(x)=Ae  '  b  )  +Ae  y  b  ) 

In  our  research  we  set  A  —  1  and  b  —  5  for  our  absorbing  boundary  conditions. 
Figure  6  depicts  the  interaction  of  the  propagating  free  particle  with  one  of  the  absorbing 
boundary  function. 
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X  X 


Figure  6.  A  propagating  wave  function  incident  on  absorbing  boundary  potential 


3.6.  Free  Particle  Propagation 

As  explained  in  Section  2.4,  the  need  to  pick  a  sufficiently  small  At  in  order  to 
avoid  complications  with  the  additive  term  in  equation  (2.13)  was  paramount.  So  in 
order  to  test  the  results  produced  from  our  numerical  propagator  we  need  to  compare 
them  with  known  analytic  solutions  of  the  same  conditions.  Since  we  have  already 
derived  the  form  of  the  free -particle  wave  function,  we  have  an  expression  with  which  to 
perform  one  of  our  first  comparisons.  Since  the  free-particle  Hamiltonian  has  no 
potential  contribution,  this  comparison  serves  to  analyze  whether  the  kinetic  component 
of  our  propagator  is  functioning  accurately.  By  checking  the  application  of  equation 
(2.19)  in  our  numerical  code  against  a  case  solely  comprising  of  one  part  of  the 
Hamiltonian  we  ensure  that  later  verification  steps  and  issues  that  may  be  found  with 
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those  verifications  are  not  attributes  of  the  kinetic  part  of  the  Hamiltonian.  Figure  7 
shows  the  comparison  between  the  numerical  result  and  the  analytic  answer  for  a  free- 
particle  propagation  derived  in  Section  2.3  and,  as  you  can  see,  they  are  in  agreement. 


t  =  5(atomic  time  units)  t  =  5(atomic  time  units) 


Figure  7.  Comparison  of  the  analytical  form  of  a  free  particle  wave  function  with  the  results  of  the  split  operator 
method.  The  blue  dots  represent  the  real  part  of  the  wave  function  as  generated  from  the  split  operator  method. 
The  lines  connecting  the  dots  represent  the  analytical  result  of  the  real  part  of  the  wave  function.  The  reason  the 
wave  function  appears  to  leave  the  window  and  not  cycle  round  again  as  discussed  in  this  section  is  that  the 
graphical  depiction  only  shows  a  narrow  portion  (-20  to  50)  of  our  x  domain. 


3.7.  Coherent  State  of  a  Harmonic  Oscillator 

After  reviewing  the  accuracy  of  our  propagator  to  a  Hamiltonian  with  no  potential 
component,  we  then  need  to  look  for  an  analytic  wave  function  for  a  Hamiltonian  with 
both  kinetic  and  potential  contributions  to  use  as  a  comparison  for  our  propagator.  The 
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wave  function  we  utilized  for  this  comparison  was  the  Coherent  State  of  the  Quantum 
Harmonic  Oscillator  [13]  given  in  the  equations  below 


/TTUUn1/4  ia)t  i(P(t ))  _fx-(x(t)) 
ip(x,  t)  =  eWa  y— — J  e  2  e  h  e  V  2A* 

(X(t))=  \—(Re[a0e-ia,t]) 

(P(t))  =  V2mhm(/?n[a0e_lwt]) 

(«*(0)2-(o:(t))2 
el0“  =  e  4 


1  /  i 

a(t)  =-ljSx(t)+-p(t) 


1  /  i 

a°  =  + 


The  expression  is  the  eigenfunction  for  the  following  Hamiltonian 


(3.9) 

(3.10) 

(3.11) 

(3.12) 

(3.13) 

(3.14) 


h2  d2  1 

H(x)  =  -  -  7  +-mco^xz 

2  max1  2 


(3.15) 


This  is  clearly  a  robust  and  complicated  expression  for  the  propagation  of  a  wave 
function  bound  within  a  quantum  harmonic  oscillator  potential.  As  such,  careful 
accounting  needs  to  be  applied  when  generating  the  analytic  result  for  our  comparison. 


3.7.1.  Ground  State 

We  first  look  at  the  simplest  case  of  the  Coherent  State,  x(0)  =  0  and  p(0)  =  0, 
as  a  comparison  for  the  propagated  wave  function  given  by  our  code.  This  is  known  as 
the  ground  state  of  the  quantum  harmonic  oscillator.  In  this  state,  the  wave  function 
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looks  very  similar  to  the  ground  state  of  the  infinite  square  well  and  its  real  and 

71 

imaginary  components  oscillate  at  a  -  phase  shift  from  each  other.  Figure  8  shows  the 

comparison  between  the  numerical  result  and  the  analytic  answer  for  the  ground  state 
coherent  function  and,  as  you  can  see,  they  are  in  agreement. 


t  =  12  t  =  18 


Figure  8.  Comparison  of  the  Coherent  function  ground  state  with  the  results  from  the  split  operator  method. 
The  blue  dots  represent  the  imaginary  part  of  the  wave  function  and  the  red  dots  represent  the  real  part  of  the 
wave  function  as  generated  from  the  split  operator  method.  The  lines  connecting  the  dots  represent  the 
analytical  result  and  the  colors  are  the  converse  for  distinction  between  the  lines  and  dots.  In  the  analytical  case 
the  red  line  represents  the  imaginary  part  and  the  blue  represents  the  real  part. 


The  comparison  of  the  free -particle  propagation  and  the  propagation  of  the 
ground  state  coherent  function  provide  a  sufficient  check  to  determine  whether  the 
numerical  propagator  of  equation  (2.19)  is  working  accurately.  Future  applications  of  the 
propagator  will  only  need  to  apply  a  sufficiently  small  At,  in  keeping  with  the  total 
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energy  of  the  system  being  analyzed,  to  provide  an  accurate  representation  of  the  time 
propagation  of  the  wave  function. 

3.8.  Discrete  Wigner  Distribution  (DWD)  Function 

Just  as  we  had  to  modify  our  approach  to  Fourier  Transforms  based  on  the 
discrete  nature  of  our  numerical  approach,  so  to  do  this  we  have  to  revisit  the  WDF  given 
in  equation  (1.1)  and  apply  it  in  a  discrete  fashion.  In  order  to  accomplish  this  we  have 
to  review  the  inner  workings  of  the  WDF.  Fundamentally,  the  WDF  is  a  Fourier 
Transform  of  a  correlated  function  in  time.  If  we  break  the  WDF  into  a  process  of 
discrete  correlation  and  a  process  of  discrete  Fourier  Transformation  we  are  half  way 
done,  given  what  we  have  already  discussed  in  Section  3.3  and  3.4. 

The  process  of  correlation  in  the  DWD  essentially  shifts  the  function,  one 
function  of  time  from  the  bottom  to  the  top  of  a  scale  in  r  space,  while  a  conjugate 
function  of  the  same  time  function  is  simultaneously  shifted  from  the  top  to  the  bottom  of 
r  space.  With  each  shift  in  r  the  product  of  the  two  functions  is  taken  and  this  builds  up 
the  correlated  function  of  time.  If  these  functions  of  time  are  localized,  that  is  to  say  they 
only  have  relatively  non-zero  values  over  a  fixed  region  of  time,  then  the  value  of  the 
correlation  will  only  have  none  trivial  values  in  the  region  of  —tmax  <  r  <  tmax  where 
tmax  is  the  last  value  of  time  for  which  the  functional  value  is  none  trivial.  This  process 
can  be  discretely  accomplished  numerically  so  that  the  only  step  left  is  to  combine  this 
process  with  the  process  of  Fourier  Transforming  the  correlated  function  in  t  space  into 
angular  frequency  space  for  all  values  of  t. 
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3.9.  Down  Sampling 

Down  Sampling  became  a  necessary  component  of  our  FORTRAN  77  code  for 
several  reasons.  In  wanting  to  explore  square  wells  with  strong  resonant  features  we 
were  forced  to  modify  the  square  well  in  one  of  two  ways;  first  being  to  increase  the 
width  and  second  being  to  increase  the  depth.  Increasing  the  width  of  the  well  would 
require  constant  movement  of  the  starting  position  of  the  wave  packet  further  to  the  left 
of  the  well  to  avoid  starting  the  computations  with  either  the  reactant  or  product  wave 
function  overlapping  the  location  of  the  well.  This  would  force  us  to  continually  readjust 
the  domain  range  of  x  to  avoid  having  the  initial  wave  function  located  in  areas  that  did 
not  have  zero  value  for  the  potential,  be  it  from  the  well  itself  or  from  the  absorbing 
boundary  conditions,  which  would  disrupt  the  initial  state.  Additionally,  this  would  also 
run  the  risk  of  diminishing  the  k  domain,  due  to  the  correlation  between  the  x  and  k 
imposed  by  the  FFT,  to  the  point  that  it  may  not  include  all  the  kinetic  energy  values 
produced  by  the  scattering. 

Altering  our  well  depth  also  imposed  complications  as  it  required  that  reduction  in 
our  time  step  to  be  able  to  satisfy  the  imposed  restriction  of  the  Split  Operator  Method. 
However,  every  time  we  decrease  the  time  step  of  our  propagator  the  length  of  the 
correlation  array  is  increased  over  a  fixed  time  interval.  That  is  to  say  is  we  subdivide  a 
length  of  time  of  10  by  intervals  of  1  we  end  of  up  with  10  data  points,  but  is  we 
subdivide  by  .1  we  end  up  with  100  data  points.  In  turn,  this  input  correlation  array 
would  increase  the  number  of  cycles  the  DWD  would  have  to  perform.  Table  3  depicts 
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this  relationship  between  the  input  array  and  the  number  of  calculations  the  DWD  would 


have  to  perform  for  a  non-down  sampled  array. 


Table  3.  Down  sampling  values 


At 

C(t)  array  length 

#  of  DWD  cycles 

Fraction  of  .09 

.5 

200 

16,040,000 

.00585 

.1 

1,000 

2,001,000,000 

.72926 

.09 

UkF 

2,743,895,583 

1 

.05 

2,000 

16,004,000,000 

5.83259 

.01 

10,000 

2,000,100,000,000 

728.92715 

The  arrays  depict  a  total  time  of  100  and  the  inclusion  of  .09  was  done  to  show  the 
attributes  for  the  value  later  utilized. 

We  decided  to  forgo  altering  the  well  width  and  instead  focused  on  changing  well 
depths,  but  to  avoid  the  calculation  issues  of  a  smaller  At,  the  correlation  array  was  down 
sampled  by  a  factor  of  9  which,  as  demonstrated  above,  made  the  total  time 
1 

approximately  —  the  actual  time  needed  to  generate  a  DWD  at  that  time  step.  As  a 

means  of  comparison,  a  DWD  with  the  .09  step  takes  about  20  minutes.  The  same  DWD 
with  a  .01  time  step  would  take  approximately  10  days13. 


12  In  this  case  the  time  length  is  99.99  as  in  our  code  that  correlation  was  computed  for  t  =  .09  to  t  =  100. 
The  instance  of  t  =  0  was  disregarded  as  by  construction  the  correlation  function  would  be  zero  at  this 
point. 

13  The  results  obtained  the  from  20  minute  trials  we  sufficient  enough  in  the  analysis  to  retain  all  susceptive 
information  of  the  WDF  that  any  increase  in  granularity  would  add  no  benefit  and  cost  excessive  time. 
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IV.  Results 


In  this  section  a  large  portion  of  the  results  obtained  from  the  code  and  methods 
described  in  Sections  II  and  III  are  displayed  for  the  benefit  of  the  reader.  Additional  results  are 
displayed  in  Appendix  A.  The  results  are  presented  in  such  a  way  as  to  give  the  reader  an 
appreciation  of  the  outcomes  of  the  modeling  and  simulations  preformed  while  saving  all 
discussion  of  analysis  and  interpretations  for  Section  V. 

4.1  Correlation  Function 

The  primary  goal  of  this  research  was  to  generate  a  Wigner  Distribution  of  the  time 
correlation  function,  given  by  equation  (1.2),  between  a  reactant  and  product  wave  function.  As 
such  it  is  necessary  to  demonstrate  the  resulting  calculated  correlation  function  obtained  via  the 
Split  Operator  Method,  already  discussed  previously,  before  moving  on  to  the  resulting  Wigner 
Distributions  obtained.  Figure  9  shows  such  a  correlation  function  generated  after  a  reactance 
wave  function  has  passed  through  a  quantum  square  well  of  depth  20Eh  and  width  of  a  =  10 a0. 


5  10  15  20  25  30  35 

t 


Figure  9.  A  correlation  function,  both  its  real  and  imaginary  components,  generated  by  the  Fortran  77  code 
utilizing  the  Split  Operator  Method.  Though  only  displayed  from  t  =  5  to  t  =  35  the  correlation  function  was 
captured  for  t  =  0.09  to  t  =  99.99 
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As  you  can  see  from  figure  9  the  sampling  interval,  as  discussed  in  the  previous  section, 
provide  more  than  sufficient  date  points  to  accurately  depict  the  function.  Additionally,  the 
correlation  function  it  given  sufficient  amount  of  time  to  reduce  to  zero  as  show  in  figure  9  which 
only  shows  a  window  of  the  function  from  t  =  5  to  t  =  35  in  order  to  show  depth  of  features  in  the 
graph.  The  true  length  of  the  correlation  function,  as  mentioned  in  the  previous  section  goes  from 
t  =  .09  to  t  =  99.99. 

4.2  Wigner  Distribution  Plots 

Using  the  method  and  code  described  in  Sections  II  and  III  respectively,  we  were 
able  to  generate  Wigner  distribution  plots  for  a  variety  of  well  depths  and  well  shapes. 
We,  of  course,  first  developed  WDF  plots  for  a  range  of  square  well  depths.  First  the 
series  between  well  depths  was  kept  small,  going  from  a  value  of  2  to  a  value  of  8  by 
increments  of  2.  After  that  point,  steps  of  20  were  used  in  order  to  reach  levels  which 
would  produce  higher  levels  of  resonant  features.  After  the  full  range  of  square  well 
DWDs  were  generated  the  same  series  of  well  depths  were  applied  to  wells  of  differing 
shapes  (triangular,  circular,  and  Gaussian).  As  means  of  a  baseline  for  our  review  of  the 
WDF  plots  generated,  figure  10  shows  the  WDF  plot  for  a  propagation  free  particle  with 
no  well  to  scatter  upon. 
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Figure  10.  The  WDF  of  a  free  particles  correlation  function 


The  unique  features  in  figure  10  are  interesting  in  and  of  themselves  and  provided  the 
overall  characteristic  seen  in  many  of  the  Wigner  Distribution  Plots  for  cases  where  a  quantum 
potential  well  was  present.  In  all  cases  the  primary  large  mountain  like  feature  seen  in  figure  10 
is  present  along  with  the  subsequent  smaller  L-shaped  ridges.  These  ridge  lines  represent  large 
portions  of  the  reactant  wave  function  correlation  with  the  product  wave  function.  Since  we  set 
out  initial  momentum  of  reactant  wave  function  relatively  large  in  order  to  avoid  spreading  the 
reactant  wave  in  the  negative  direction  we  gain  the  effect  of  supplying  a  large  amount  of  initial 
energy  into  the  correlation  and  thus  in  all  of  our  Wigner  Distributions  this  very  large  primary 
peak  is  present. 

We  have  included  the  WDF  plots  of  the  square  wells  and  circular  wells  in  this 
section.  Appendix  A  displays  the  results  of  the  triangular  and  Gaussian  wells  due  to  their 
lack  of  structure.  A  full  analysis  of  these  results  is  found  in  Section  V. 
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As  described  in  Section  2.5  of  the  square  well  scattering,  the  potential  for  a  square 


well  is  a  set  of  functions  as  show  in  equation  (4.1) 


—V0>  —a<x<a 
\x\  >  a 


™  -  n: 

This  is  the  form  used  for  the  Vint  term  of  equation  (3.3). 


(4.1) 


Since  the  potential  is  discretely  defined,  the  potential  is  an  approximation  of  the 
square  well.  The  smaller  we  make  Ax,  the  closer  we  make  the  approximation  to  a  true 
square  well,  but  ultimately  the  numerical  process  can  never  truly  replicate  the  piecewise 
function  described  in  equation  (4.1).  Fortunately  the  distortion  is  not  large  enough  to 
cause  issue  in  our  data  with  the  Ax  we  have  selected,  Ax  =  .0976. 
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4.2.1.  Square  Well  with  depth  V„  =  2,  4,  6,  &  8  Eh 


Figure  11.  Depiction  of  a  WDF  plot  generated  from 
the  time  correlation  function  of  a  reactant  wave 
function  scattering  through  a  quantum  potential  well 
of  2Eh.  The  red  lines  indicate  frequencies  of  total 
transmission  that  correspond  to  a  well  of  that  depth. 


Figure  12.  Depiction  of  a  WDF  plot  generated  from 
the  time  correlation  function  of  a  reactant  wave 
function  scattering  through  a  quantum  potential  well 
of  4Efl.  The  red  lines  indicate  frequencies  of  total 
transmission  that  correspond  to  a  well  of  that  depth. 


Figure  13.  Depiction  of  a  WDF  plot  generated  from 
the  time  correlation  function  of  a  reactant  wave 
function  scattering  through  a  quantum  potential  well 
of  6 Eh.  The  red  lines  indicate  frequencies  of  total 
transmission  that  correspond  to  a  well  of  that  depth. 


Figure  14.  Depiction  of  a  WDF  plot  generated  from 
the  time  correlation  function  of  a  reactant  wave 
function  scattering  through  a  quantum  potential  well 
of  8 Eh.  The  red  lines  indicate  frequencies  of  total 
transmission  that  correspond  to  a  well  of  that  depth. 


42 


4.2.2.  Square  and  Circular  Well  with  depth  Vo  -10,  20,  40,  60,  &  80  Eh 

The  FORTRAN  77  code  generated  was  tailored  to  the  well  width,  2a,  of  20 
(atomic  units)  for  the  square  well.  When  the  code  was  modified  to  examine  the  circular 
well  the  following  equation  was  used 


(4.2) 


V(x) 


)  —  \V02 - x2,  —a<x<a 


a 


0,  |x|  >  a 

This  is  the  form  used  for  the  Vint  term  of  equation  (35).  In  this  form  the  circular  well 
maintains  the  same  overall  spacing  in  x  as  the  square  well  did,  so  that  the  new  potential 
can  simple  be  substituted  in  the  pre-existing  FORTRAN  77  code. 


Figure  15.  Depiction  of  a  WDF  plot  generated  from 
the  time  correlation  function  of  a  reactant  wave 
function  scattering  through  a  quantum  square  well  of 
10 Eh.  The  red  lines  indicate  frequencies  of  total 
transmission  that  correspond  to  a  well  of  that  depth. 


Figure  16.  Depiction  of  a  WDF  plot  generated  from 
the  time  correlation  function  of  a  reactant  wave 
function  scattering  through  a  quantum  circular  well  of 
10  Eh. 
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Figure  17.  Depiction  of  a  WDF  plot  generated  from 
the  time  correlation  function  of  a  reactant  wave 
function  scattering  through  a  quantum  square  well  of 
20  Eh.  The  red  lines  indicate  frequencies  of  total 
transmission  that  correspond  to  a  well  of  that  depth. 


Figure  18.  Depiction  of  a  WDF  plot  generated  from 
the  time  correlation  function  of  a  reactant  wave 
function  scattering  through  a  quantum  circular  well  of 
20  Eh. 


Figure  19.  Depiction  of  a  WDF  plot  generated  from 
the  time  correlation  function  of  a  reactant  wave 
function  scattering  through  a  quantum  square  well  of 
40  Eh.  The  red  lines  indicate  frequencies  of  total 
transmission  that  correspond  to  a  well  of  that  depth. 


Figure  20.  Depiction  of  a  WDF  plot  generated  from 
the  time  correlation  function  of  a  reactant  wave 
function  scattering  through  a  quantum  circular  well  of 
40  Eh. 
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Figure  21.  Depiction  of  a  WDF  plot  generated  from 
the  time  correlation  function  of  a  reactant  wave 
function  scattering  through  a  quantum  square  well  of 
60  Eh .  The  red  lines  indicate  frequencies  of  total 
transmission  that  correspond  to  a  well  of  that  depth. 


Figure  22.  Depiction  of  a  WDF  plot  generated  from 
the  time  correlation  function  of  a  reactant  wave 
function  scattering  through  a  quantum  circular  well  of 
60  Eh. 


Figure  23.  Depiction  of  a  WDF  plot  generated  from 
the  time  correlation  function  of  a  reactant  wave 
function  scattering  through  a  quantum  square  well  of 
80  Eh.  The  red  lines  indicate  frequencies  of  total 
transmission  that  correspond  to  a  well  of  that  depth. 


Figure  24.  Depiction  of  a  WDF  plot  generated  from 
the  time  correlation  function  of  a  reactant  wave 
function  scattering  through  a  quantum  circular  well  of 
80  Eh. 
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V.  Analysis 


In  viewing  the  results  of  the  square  well  and  comparing  those  results  with  known 
attributes  of  that  same  well,  many  noteworthy  features  begin  to  appear  in  our  WDF  plots. 
We  will  examine  these  attributes  in  the  case  of  a  square  well  with  a  depth  of  20 Eh  unless 
otherwise  noted  in  the  discussion.  This  is  not  to  say  that  these  features  are  not  present  in 
all  cases,  but  the  example  gives  us  a  good  case  where  we  do  not  run  into  extremes  of  all 
the  features,  which  can  obscure  the  individual  attributes.14 

5.1.  Transmission  Coefficients 

As  discussed  in  Section  2.5.1  a  square  well  allows  for  certain  energy  bands  to 
pass  through  as  if  the  well  were  transparent.  When  we  take  the  integrate  projection  of  the 
WDF  generated  from  the  scattering  event  towards  the  frequency  axis,  we  get  the 
probability  distribution  of  the  frequency  of  our  correlation  function  from  that  scattering 
event.  Figure  25  shows  the  integrate  projection  of  the  WDF  for  the  square  well  of 
V0  —  20  and  a  —  10. 


14  This  also  helps  us  from  diluting  the  results  with  multiple  graphical  plots  restating  the  point. 
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Figure  25.  The  integrated  projection  of  the  WDF  in  frequency  is  given  by  the  blue  dots  and  the  transmission 
function  of  a  square  well  of  width  2a  =  10  and  V0  =  20  is  given  by  the  solid  blue  line. 


By  inspection,  the  distribution  has  an  overall  Gaussian  shape.  This  is  due  to  the  fact  that 
the  reactant  and  product  wave  functions  were  Gaussian  by  construction15.  Within  the 
overall  Gaussian  shape  there  appears  to  be  protruding  peaks  at  specific  values  of 
frequency.  These  peaks  coincide  when  compared  to  the  known  transmission  peaks  of  the 
square  well.  Further,  the  overall  Gaussian  function  is  given  by  the  modulus  of  the 
frequency  correlation  of  the  two  initial  wave  functions  when  no  well  is  present 

|  C  (eo)  | 2  =  C*(oo)C((o)  (5.1) 


15  If  we  had  use  Lorentzian  function  then  the  overall  shape  of  the  correlation  would  be  Lorentzian. 
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Which  in  our  particular  case  of  x0  =  —20 ,k0  =  3  and  8  =  1  is  the  following 


|C(<d)I 


le-4(9+20V-2(U)+2cj)  ^  _|_  e(80-24i)V-2w  _j_  g(80+24i)V-2cJ  _)_  gl60V2oj^  ^  . 


The  production  of  equation  (5.2)16  with  the  transmission  function  (2.27)  of  the  square 
well  gives  us  the  same  result  as  our  integrated  projection  of  the  WDF  upon  the  frequency 
axes.  Figure  26  shows  a  comparison  of  the  results  of  the  WDF  integrated  projection  upon 
the  frequency  axis  and  the  product  of  the  transmission  function  with  equation  (5.2) 


Figure  26.  The  comparison  of  the  product  of  the  transmission  function  (2.27)  and  equation  (5.2)  and  the  result  of 
the  integrated  projection  of  the  WDF  upon  the  frequency  axis.  The  red  line  is  the  product  result  and  the  blue 
dots  are  the  results  of  the  integration.17 


16  This  result  was  obtained  by  taking  the  Fourier  Transform  of  C  (t)  for  the  case  where  there  was  no 
potential  well  involved  and  finding  the  modulus  of  the  result  in  Mathematica  6.0. 

17  The  slight  differences  in  the  red  line  and  the  blue  dotted  line  potential  arise  from  the  discrete  nature  of 
the  potential  well  and  the  fact  that  as  mentioned  in  Section  4.2  the  well  is  actually  trapezoidal.  It  is  the 
opinion  of  the  author  however  that  the  difference  has  more  to  do  with  the  analytical  approximation 
Mathematica  6.0  made  to  reach  the  result  in  equation  (5.2). 
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The  integrated  projection  in  time  of  a  WDF  generated  from  the  time  correlation  of 
a  wave  function  scattering  off  a  potential  well  is  the  product  of  the  transmission  function 
of  that  potential  well  and  the  modulus  of  the  frequency  correlation  of  the  wave  function 
in  free  space.  This  property  can  be  expanded  to  wells  of  unknown  transmission  functions 
and  from  it  we  can  derive  the  transmission  functions  of  those  wells.  By  this  we  mean  that 
we  can  take  the  correlation  function  of  a  scatter  event  off  of  those  wells,  generate  a 
Wigner  Distribution  of  that  scattering,  and  integrate  that  distribution  in  time  to  get  a 
probability  density  in  terms  of  frequency.  From  here  we  can  divide  the  probability 
density  function  by  the  modulus  of  the  correlation  function  in  frequency  of  a  non¬ 
scattering  event  and  from  that  arrive  at  out  transmission  function,  which  is  what  we 
should  expect. 

5.2.  Square  Well  Energy  Levels 

Previously  we  saw  that  the  integrated  projection  of  the  WDF  corresponded  to  the 
transmission  function  of  the  square  well  we  took  our  time  correlation  from.  It  is  no 
surprise  then  to  see  that  the  resonant  frequency  levels  match  with  the  peak  values  of  the 
WDF’s  integrated  projection  in  frequency.  However,  if  we  take  these  frequency  lines 
and  place  them  across  our  3  dimensional  WDF  surface,  we  see  that  the  resonant 
frequency  levels  connect  the  positive  peaks  of  every  ridge  structure  in  our  surface.  This 
is  shown  in  figure  27. 


49 


W(U>) 


10 


Figure  27.  WDF  plot  with  total  transmission  energy  level  bands  of  a  well  width  of  40Fft.  The  red  lines  indicate 
allowed  frequencies  for  a  quantum  square  well  of  depth  40 Eh  and  width  of  20a0. 


However,  that  is  not  to  say  that  these  lines  correspond  to  every  positive  peak 
position.  If  we  examine  the  peak  values  of  the  third  ridge  structure  in  figure  27,  by  third 
ridge  we  mean  the  third  positive  structure  seen  as  you  progress  along  the  time  axis,  we 
see  that  the  frequency  levels  correspond  to  every  other  peak  structure  in  that  ridge  line. 
There  appears  to  be  peaks  that  correspond  to  half  way  frequencies  of  the  resonant 
frequencies.  When  the  frequency  levels  of  a  well  with  double  the  actual  scattered  wells 
width  are  overlaid  on  top  of  the  WDF  surface  we  see  that  all  the  peaks  in  the  third  ridge 
line  correspond  to  the  new  resonant  frequencies,  but  we  also  notice  that  the  preceding 
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ridge  line,  or  second  as  you  progress  along  the  time  axis,  contains  negative  valley 
position  along  these  lines  “half’  frequencies.  This  is  shown  in  figure  28. 
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Figure  28.  WDF  plot  with  total  transmission  energy  level  bands  of  a  well  width  of  40 Eh.  The  red  lines  indicate 
allowed  frequencies  for  a  quantum  square  well  of  depth  40 Eh  and  width  of  40a0. 


Clearly  the  summation  of  these  “half’  frequencies  is  what  leads  to  the  positions  of 
anti-resonant  in  the  WDF’s  integrate  projection  in  frequency,  but  the  meaning  behind 
these  negative  values  and  corresponding  positive  features  is  unclear. 

5.3.  Time  Arrivals 

Yet  another  interesting  feature  we  can  observe  in  the  WDF  plot  is  the  time  arrival 
of  a  classical  particle  with  an  energy  defined  by  h<y  in  comparison  with  what  we  are 
seeing  in  the  WDF  plot.  In  the  case  of  a  classical  particle  with  initial  position  equal  to 
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the  initial  position  of  the  wave  function  and  energy  equal  to  h a>,  we  can  determine  its 
arrival  to  the  left  edge  of  the  square  well  by  the  following  expression 


Where  x  —  b  and  v  are  defined  by  the  following, 


(5.3) 


£"/,///  =  h  co  =  -mv2 


Since  the  time  it  would  take  to  travel  to  and  travel  from  the  well  is  the  same,  we  have 
denoted  equations  (5.4)  and  (5.5)  for  both  region  I,  the  region  before  the  well,  and  region 
III,  the  region  after  the  well.  Once  the  particle  enters  the  well,  we  assume  it  moves  from 
the  left  edge  to  the  right  at  a  velocity  described  by  the  following  expressions 


1  _ 

Eu  =  h  co  —  -  mv 2  —  V0 

vn  = 

Within  the  well  the  particle  can  do  many  different  things.  It  can  pass  through  the  well.  It 
can  bounce  off  the  right  edge  of  the  well  and  exit  the  well  through  the  left  edge.  It  can 
bounce  off  the  right  edge  and  then  bounce  off  the  left  edge  and  exit  the  well  back  through 
the  right  edge.  As  you  can  see,  the  particle  has  an  integer,  n,  number  of  scenarios  of 
travel  within  the  well  before  it  finally  exits.  These  scenarios  can  be  expressed  by 
multiplying  the  width  of  the  well  by  the  number  of  time  the  particle  traverses  the  well 
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before  exiting,  n.  A  summation  of  all  these  expressions  is  given  in  equation  (5.8)  where 
the  total  time  it  takes  the  particle  to  reach  the  position  of  the  stationary  wave  function  is 
determined. 


t(m)  = 


2b 

; - +  71 


2a 


(5.8) 


I2ha>  12  (h  a>  +  V0 ) 
a/  m  \  m 

If  we  plot  equation  (5.8)  for  each  positive  integer  value  of  n  and  set  the  resulting 
lines  over  top  of  our  WDF  plots,  as  done  in  figure  29,  we  can  see  that  the  overall  ridge 
structure  seen  in  the  WDF  plots  correspond  to  this  classical  analogy. 
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Figure  29.  A  contour  plot  of  the  WDF  with  yellow  lines  indicating  the  arrival  times  of  a  classical  particle  with 
the  given  frequency/energy  while  the  red  line  indicate  allowed  frequencies  of  the  quantum  well. 
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While  the  first  line  represents  the  case  where  the  classical  particle  proceeded 
directly  through  the  square  well,  the  second  line  represents  the  case  where  the  particle 
bounced  once  off  of  the  right  wall  of  the  well  and  exited  the  left  side  of  the  well,  which  is 
the  side  the  particle  entered  the  well  initially.  It  is  surprising  to  see  that  it  is  this  classical 
example  which  corresponds  to  the  ridge  line  where  the  positive  peak  values  only  occur  at 
the  energy  levels  of  a  square  well  of  this  type,  (i.e.  a  =  10  and  V0  =  20  ).  Classically, 
that  information  is  leaving  the  well  in  left  direction,  but  quantum  mechanically  that 
information  is  being  recorded  in  the  right  region  of  the  well. 

5.4.  Spectrogram  versus  WDF 

Often  in  signal  processing,  the  technique  of  taking  a  spectrogram  is  used  in  order 
to  gauge  contributions  of  frequency  components  within  a  signal.  This  is  predominately 
used  to  weed  out  ambient  noise  and  bolster  desired  frequencies.  A  spectrogram  is  the 
product  of  squaring  a  Short  Time  Fourier  Transform,  depicted  in  equation  (5.9) 

Z(T.»)  =  fs(tMt-r)e-*‘#  15  9) 

Spectrogram  —  \ x(r,  co)  | 2  (5.10) 

The  windowing  function,  w(t  —  r),  can  be  a  multitude  of  functions  from  a  Hanning 
function,  a  step  function,  to  a  Gaussian. 

We  applied  this  windowed  Fourier  Transform  to  our  correlation  function,  as 
shown  in  figure  30,  in  order  to  see  what  information  it  provided  of  the  scattering  event  in 
comparison  to  what  we  have  gained  from  our  WDF  plots. 
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Figure  30.  The  windowed  FT  of  our  correlation  function  utilizing  various  FWFUVI  values  for  the  Gaussian 
window  term.  From  right  to  left  the  value  of  c  in  equation  (5.11)  is  as  follows;  4.24661,  2.1233,  1.06165, 
0.707768,  and  0.530826  respectively. 


We  used  a  Gaussian  for  our  windowing  function 

(t-T)2  (5.11) 

w(t  —  T)  =  e  2At2 

From  a  cursory  observation  of  the  plots  shown  in  figure  30  we  see  that  most  of  the 
information  returned  in  the  spectrogram  is  obscured  by  the  sure  volume  of  energy 
recorded  in  the  initial  large  bulge  of  the  correlation  function.  However,  by  applying  a 
logarithmic  scale  to  the  spectrogram  we  were  able  to  bring  detailed  elements  of  the 
scattering  event  back  up  to  a  scale  where  we  could  compare  them  with  our  WDF.  Figure 
31  shows  these  revised  plots  of  the  spectrograms. 
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Figure  31.  A  series  of  spectrogram  plots  of  a  correlation  function  of  the  scattering  event  involving  a  well  of 
depth  40 Eh  and  width  of  20a0.  The  box  indicates  the  size  in  comparison  to  the  overall  time  domain  that  the 
Gaussian  window  function,  equation  (5.11),  set  to  when  the  spectrogram  was  generated.  A  Nth  window  size 
corresponds  to  a  At  given  by  the  overall  time  of  the  correlation  function  divided  by  N,  AtN  =  99. 99/N. 


57 


As  shown  in  figure  31,  a  series  of  spectrograms  were  taken  for  the  case  of  a 
correlation  function  generated  from  the  scattering  through  a  quantum  potential  square 
well  of  depth  40 Eh  and  width  of  20<x0 .  Each  new  spectrogram  had  the  window  size  of  the 
gaussian  window  function,  equation  (5.1 1),  reduced  in  relation  to  the  overall  time  of  the 
correlation  function.  Since  the  length  of  time  of  the  correlation  function  was  99.99(atomic  time 
untis)  the  gaussian  window  function  was  set  to  first  be  a  20th  of  that  size,  4.9995(atomic  time 
units),  then  a  30th,  3.333(atomic  time  units),  and  so  on  until  a  window  size  that  was  a  100th  the 
time  scale,  .999  atomic  time  units.  As  was  done  with  the  WDF  plots,  the  allowed  frequencies  and 
time  arrival  information  was  enterposed  over  top  the  contour  plots  of  the  spectrograms. 

As  you  can  see  the  allowed  frequency  information  can  be  determined  in  spectrograms 
where  a  wide  window  was  taken,  20th  through  about  40th.  However,  as  the  windowing  function  is 
narrowed,  or  we  gain  a  greater  degree  of  clarity  in  time,  we  loose  our  information  in  frequency. 

In  figure  31,  the  40th  window  spectrogram  starts  to  show  that  the  null  regions,  the  small  blue 
circles,  between  the  allowed  frequency  bands  begin  to  drift  away  from  the  center  of  those 
frequency  bands  until,  as  seen  in  the  50th  window  spectrogram,  the  null  is  located  in  regions  of 
allowed  frequency  and  in  the  60th  window  spectrogram  the  null  has  moved  into  another  spacing 
between  energy  bands. 

On  the  converse  as  we  loose  this  clarity  frequency  in  smaller  windowed  spectrograms  we 
tend  to  gain,  to  a  small  extent,  clarity  in  time.  All  the  spectrograms  have  the  same  general  sweep 
as  our  time  arrival  lines,  but  it  is  not  until  you  reach  the  40th  window  spectrogram  and  then  again 
in  the  50th  window  spectrogram  where  you  see  the  shape  and  place  meant  of  these  lines  in  nearest 
agreement.  However,  this  is  precisely  the  window  scales  that  positioned  out  null,  anti-resonance, 
at  a  position  of  resonance. 
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The  spectrogram  is  a  helpful  tool  in  signal  processing  since  the  spectrogram  is  always  a 
positive-valued  distribution  of  energy  in  the  signal.  However  the  spectrogram  comes  at  a  cost 
that  requires  you  to  make  a  tradeoff  in  the  precision  of  your  time  or  frequency  information.  If  we 
compare  the  outcomes  of  the  spectrogram  with  our  WDF  shown  in  figure  29  we  find  that  the 
WDF  enumerates  the  same  information  (as  shown  in  Section  5.2  and  Section  5.3),  but  requires  no 
trade-off  be  made  in  the  time  or  frequency18. 

5.5.  Conclusion 

The  application  of  a  correlation  function  to  a  scattering  event  and  the  production 
of  its  WDF  have  yielded  a  wealth  of  information  about  the  long  term  effects  a  potential 
well  has  on  a  wave  function.  By  analyzing  the  WDF  plots  of  scattering  events  we  have 
seen  that  is  it  possible  to  locate  the  allowed  frequencies  of  the  potential  well. 
Additionally,  the  WDF  tells  us  something,  if  only  in  terms  of  a  classical  kinematic 
analogy,  of  the  time  arrival  of  information  of  the  reactant  wave  function  after  it  has 
passed  through  the  well. 

Furthermore,  the  process  of  taking  WDF  plots  of  the  correlation  function  of  a 
scattered  reactant  wave  function  gives  us  a  set  of  tools  to  analyze  the  reaction  of  quantum 
bodies  with  almost  any  potential  well.  Some  of  these  tools  include  the  reconstruction  of 
the  particular  potential  well’s  transmission  function  as  well  as  the  time  departure  beyond 
the  potential  frequencies.  While  the  negative  attributes  of  the  WDF  remain  unintuitive  to 
analysis,  the  WDF  provides  more  insight  into  the  potential  scattering  event  then  is  gained 
from  a  Spectrogram  of  the  same  event.  In  the  cases  demonstrated  here,  the  WDF  can  be 

18  Although  the  ability  for  the  WDF  to  hold  negative  values  encourages  us  to  avoid  looking  at  it  as  an 
energy  distribution  as  we  were  able  to  in  the  spectrogram,  it  is  clearly  analogous  in  most  respects. 
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used  in  place  of  a  spectrogram  to  return  the  same  information  the  spectrogram  could 
provide  without  the  need  for  scaling  window  size.  It  is  in  the  opinion  of  the  author  that 
applications  of  WDF  plots  in  broader  chemical  reactions  will  provide  researchers  in  the 
field  of  chemical  physics  a  broader  interpretation  of  the  reaction  as  a  whole. 

5.6.  Further  Research 

Recommendations  of  further  work  would  be  to  generate  correlation  functions  of 
the  exiting  propagating  wave  function  from  both  sides  of  the  potential  well  and  compare 
their  WDF  plots.  As  well  as  comparing  the  WDF  plots  to  other  signals  processing 
techniques  such  as  the  ambiguity  function  (see  Appendix  B).  Additionally,  expanding 
the  application  of  generating  WDF  from  correlation  functions  of  2  and  3  dimensions 
potential  scattering  events  as  well  as  coupled  potential  wells. 
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Appendix  A 


The  WDF  produced  by  scattering  our  propagating  wave  through  a  series  of  Gaussian  and 
triangular  potential  well  showed  remarkably  little  results.  Since  the  FORTRAN  code  was 
designed  to  have  the  potential  well  occupy  the  region  of  x  =  -10  to  x  =  10  it  is  possible  that  as  we 
increased  the  depth  of  the  potential  well,  the  overall  shape  of  our  Gaussian  function  and  triangular 
functions  became  closer  approximations  of  a  delta  function  then  occurred  in  the  square  or  circular 
case.  In  the  case  of  a  delta  function  there  is  no  internal  well  width  for  the  wave  to  bounce  within 
and  the  propagation  become  analogous  to  the  free  particle  example  with  no  well  to  scatter  upon. 
When  the  wells  are  shallower  then  it  is  also  possible  that  the  /?  term  is  not  large  enough  for 
resonant  feature  to  become  prominent  over  the  initial  transfer  of  frequency  as  the  wave  pass 
through  the  wall  once. 

Additionally,  it  is  possible  that  the  angle  of  the  well  barrier  is  what  cased  this  lack  of 
feature.  In  the  Square  and  Circular  well  the  angle  made  between  the  region  of  space  with  no  well 
and  the  region  of  space  with  a  well  is  90°.  In  the  case  of  the  Gaussian  and  Triangular  Wells  that 
angle  is  either  gradual  or  less  than  90°.  This  is  depicted  in  figure  3 1 . 


U  \J 

Figure  32.  The  well  shapes  used  as  the  potential  component  of  the  Spilt  Operator  Method.  From  right  to  left 
they  are  a  square  well,  a  circular  well,  a  triangular  well,  and  a  Gaussian  well. 
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Gaussian  Well 


Provided  below  are  the  WDF  plot  generated  for  a  potential  well  given  by  equation 
(A.l)  where  V0  is  the  depth  of  the  well. 


V(x)  =  -V0e~x2 


(A.l) 


Gaussian  V0  =  10  Eh 


Figure  33.  Depiction  of  a  WDF  plot  generated  from  the  time  correlation  function  of  a  reactant  wave  function 
scattering  through  a  quantum  Gaussian  well  of  10 Eh. 

Gaussian  V0  =  20  Eh 


Figure  34.  Depiction  of  a  WDF  plot  generated  from  the  time  correlation  function  of  a  reactant  wave  function 
scattering  through  a  quantum  Gaussian  well  of  20 Eh. 
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Gaussian  V0  =  40  Eh 


Figure  35.  Depiction  of  a  WDF  plot  generated  from  the  time  correlation  function  of  a  reactant  wave  function 
scattering  through  a  quantum  Gaussian  well  of  40 Eh. 


Gaussian  V0  =  60  Eh 


Figure  36.  Depiction  of  a  WDF  plot  generated  from  the  time  correlation  function  of  a  reactant  wave  function 
scattering  through  a  quantum  Gaussian  well  of  60 Eh. 
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Gaussian  V0  =  80  Eh 


Figure  37.  Depiction  of  a  WDF  plot  generated  from  the  time  correlation  function  of  a  reactant  wave  function 
scattering  through  a  quantum  Gaussian  well  of  80 Eh. 


Triangular  Well 

Provided  below  are  the  WDF  plot  generated  for  a  potential  well  given  by  equation 
(A. 2)  where  V0  is  the  depth  of  the  well  and  a  is  the  width  of  the  well  at  the  base  of  the 
triangle. 


V(x) 


Vb 

l-Vo  +  ^lxl 

0, 


—a  <  x  <  a 
x\  >  a 


(A.2) 
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Triangular  Va  =  10  Eh 


Figure  38.  Depiction  of  a  WDF  plot  generated  from  the  time  correlation  function  of  a  reactant  wave  function 
scattering  through  a  quantum  triangular  well  of  10 Eh. 


Triangular  V0  =  20  Eh 


Figure  39.  Depiction  of  a  WDF  plot  generated  from  the  time  correlation  function  of  a  reactant  wave  function 
scattering  through  a  quantum  triangular  well  of  20 Eh. 
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Appendix  B 


In  addition  to  rescaling  the  spectrogram  data  by  use  of  taking  the  logarithmic 
value  at  each  point  in  the  time  frequency  domain  we  also  normalized  those  same  values 
versus  the  entirety  of  the  spectrogram.  Be  taking  this  approach  we  were  able  to  bring  up 
smaller  features  that  were  obscured  by  the  size  of  the  primary  peaks  in  the  spectrogram. 
Below  (Figure  40  and  41)  is  the  results  of  these  normalized  spectrogram  displayed  in  a 
similar  manner  to  the  plots  shown  in  figure  3 1 . 


Figure  40.  Normalized  spectrograms  of  the  generated  correlation  function  of  the  scattering  event  involving  a 
well  of  depth  40 Eh  and  width  of  20ao.  The  red  lines  indicate  the  classical  time  arrival  of  a  particle  with  energy 
given  by  <d  and  the  yellow  lines  represent  the  allowed  transmission  resonances  of  the  square  well.  Window  sizes 
of  a  10th,  20th,  30th  and  40th  of  the  overall  time  length  of  the  correlation  are  given  shown. 
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Figure  41.  Normalized  spectrograms  of  the  generated  correlation  function  of  the  scattering  event  involving  a 
well  of  depth  40 Eh  and  width  of  20 a0.  The  red  lines  indicate  the  classical  time  arrival  of  a  particle  with  energy 
given  by  a  and  the  yellow  lines  represent  the  allowed  transmission  resonances  of  the  square  well.  Window  sizes 
of  a  50lh,  60th,  70lh  and  80,h  of  the  overall  time  length  of  the  correlation  are  given  shown. 


We  observed  that  once  the  normalized  spectrogram  was  taken  it  was  far  easier  to 
draw  a  comparison  between  the  classical  time  arrival  analysis  and  the  information  seen  in 
the  spectrogram.  However,  the  tradeoff  between  the  accuracy  of  the  frequency 
component  to  the  time  component  of  the  spectrogram  is  still  present  as  well  saw  in  the 
logarithmic  version.  As  we  refine  our  windowing  function  in  time  we  can  see  that  the 
null  regions  between  out  transmission  resonance  lines  move  out  of  the  regions  between 
the  resonance  lines  and  spread  to  overlap  regions  of  transmission.  As  discussed 
previously,  this  time/frequency  trade-off  is  absent  in  the  Wigner  Distribution  Function. 
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Appendix  C 


Another  function  used  in  signal  processing  is  the  Ambiguity  function  as  shown  in 
equation  (B.l).  [19] 

x(t,  <x>)  =  J  s(t  +  r)s*(t  —  r)e_lftdt 

This  function  is  typically  used  in  the  radar  community  to  observe  the  Doppler  delay  of  a  radar 
pulse  upon  a  moving  target.  Figure  40  depicts  the  real  and  imaginary  components  of  the 
Ambiguity  function  for  a  complex  Gaussian  centered  about  x  =  2.5  and  all  other  attributes  set  to 
unity.  Figure  41  depicts  the  modulus  of  the  ambiguity  function  of  the  same  Gaussian  already 
described. 
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Figure  42.  Ambiguity  function  (real  component  on  the  left  and  imaginary  component  on  the  right)  of  a  complex 
Gaussian  centered  about  x  =  2.5  and  all  other  attributes  set  to  unity. 
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Figure  43.  The  modulus  of  a  Ambiguity  function  of  a  complex  Gaussian  centered  about  x  =  2.5  and  all  other 
attributes  set  to  unity. 


This  function  could  also  be  used  to  analyze  the  correlation  function  generated  in  our  methods 
described  in  Section  II.  An  example  of  such  an  ambiguity  function  if  shown  in  figure  42  for  the 
case  of  a  correlation  function  generated  from  scattering  through  a  quantum  potential  square  well 
of  depth  40 Eh  and  width  of  20a0. 
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Figure  44.  The  modulus  of  an  Ambiguity  Function  of  a  correlation  function  generated  from  scattering  through  a 
quantum  potential  square  well  of  depth  40Eh  and  width  of  20a0.  Represented  on  a  lograthmic  scale. 
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